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Abstract 

The turbulent mixing noise of a supersonic jet is 
calculated for a round convergent-divergent nozzle at 
the design pressure ratio. Aerodynamic computations 
are performed using the PARC code with a k-c tur- 
bulence model. Lighthill’s acoustic analogy combined 
with Ribner’s assumption is adopted. The acoustics so- 
lution is based upon the methodology followed by GE 
in the MGB code. The source correlation function is 
expressed as a linear combination of second-order ten- 
sors. Assuming separable second-order correlations and 
incorporating Batchelor’s isotropic turbulence model, 
the source term was calculated from the kinetic en- 
ergy of turbulence. A Gaussian distribution for the 
time-delay of correlation was introduced. The CFD so- 
lution was used to obtain the source strength as well 
as the characteristic time-delay of correlation. The ef- 
fect of sound/flow interaction was incorporated using 
the high frequency asymptotic solution to Lilley’s equa- 
tion for axisymmetric geometries. Acoustic results in- 
clude sound pressure level directivity and spectra at 
different polar angles. The aerodynamic and acoustic 
results demonstrate favorable agreement with experi- 
mental data. 

Introduction 

Recent advances in the area of computational fluid 
dynamics have provided researchers with new tools for 
approaching aeroacoustic problems. Though the di- 
rect computation of the sound field and its propaga- 
tion mechanism is the ultimate goal in an area which 
has been labeled as computational aeroacoustics, the 
extraordinary difficult nature of these computations at 
the present time lead to a more practical approach, 
namely using the available computational resources. 
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Such an approach is taken here to study the high speed 
jet noise because of the current interest in High Speed 
Civil Transport (HSCT). 

Over a decade ago, a unified aerodynamic/acoustic 
prediction model was developed by GE for the Federal 
Aviation Administration [2]. The aerodynamic pre- 
dictions within this code are carried out by applying 
an extension of Reichardt’s model. The three compo- 
nents of turbulent shearing stress are computed based 
upon contour integrals around the nozzle exit geome- 
try and utilized in deriving an expression for the source 
strength. Though the closed form nature of Reichardt’s 
solution has led to a versatile and relatively fast numer- 
ical scheme, the assumptions inherent in the derivations 
render inaccurate aerodynamic prediction for complex 
geometries. 

Reichardt’s solution essentially assumes exhaust 
nozzle elements discharge axially. It neglects radial 
mean flow and swirl and does not take into account 
the effect of shock structure on mixing and turbu- 
lence. These assumptions, among others, could vitiate 
the aerodynamic predictions substantially. In the light 
of these assumptions the model proves insufficient for 
noise prediction of new lobed hypermix nozzles that 
are the present focus in HSCT research. To circumvent 
some of the pitfalls associated with Reichardt’s model, 
a two-stage algorithm is considered in which the aero- 
dynamic calculations are carried out independently and 
the resulting plume data are introduced for noise com- 
putation. The following computational advantages are 
associated with a two-stage algorithm: 

(1) the aerodynamic calculations can start from within 
the nozzle, giving a more realistic model of the exit 
conditions; 

(2) grid selection can be done independently for aero- 
dynamic and acoustic purposes. The plume is 
treated as an integration of axial slices in acous- 
tic computations; 

(3) proper prediction of shock-cell structure and inten- 
sity for off-design conditions improves shock noise 
assessment. 


1 


The geometry considered in this paper is a round 
convergent-divergent nozzle. This particular nozzle has 
previously been tested aerodynamically and acousti- 
cally for NASA-Lewis by General Electric [1]. Test data 
has been especially helpful in assessing the CFD results 
before proceeding with noise computations. The flow- 
path is designed to obtain an isentropic, uniform and 
parallel flow at the nozzle exit for the design Mach num- 
ber of 1.4, thereby eliminating to a large degree, any 
shock induced noise. An axisymmetric version of PARC 
code with Chien’s k-e turbulence model [3] has been 
used. The complete Reynolds-averaged Navier Stokes 
equations are solved in cons erva t ive law form. The vis- 
cus coefficients are determined based on Sutherland’s 
viscosity law, Stokes hypothesis and constant Prandtl 
number assumption. The jet plume flow field is com- 
pared with data as well as solution from Reichardt’s 
model. 

Method of Solution 

The solution technique is essentially based upon 
the methodology followed by GE in developing the 
MGB code [2]. Reich ardt’s aerodynamic model is re- 
placed with CFD and the technique is briefly outlined 
here with emphasis on the new changes and compar- 
isons with the original model. 

The solution technique is described in two sec- 
tions, source spectrum model and sound/flow inter- 
action. The first section describes the application of 
PARC with a k-e turbulence model to computation 
of source strength and its spectrum and explains the 
empirical constants used in computation of supersonic 
convection factor. In the second section the effect of 
surrounding medium on radiated noise from convecting 
quadrupole sources is described. It is explained how the 
mean flow temperature and velocity, as obtained from 
CFD, are used in computation of sound/flow interac- 
tion of axisymmetric jets. 


corresponding spectrum can be expressed in terms of 
the Fourier transform of the autocorrelation function 

p) 

As such, for a quasi-incompressible turbulence, the 
source strength is characterized by a two-point time- 
delayed fourth-order velocity correlation tensor. For 
a nearly parallel mean flow, the constituent velocities 
may be written as the sum of the local mean veloc- 
ity and the turbulent velocity, V* = USn + v*. Sub- 
stituting Vi into equation (1), the contribution to the 
integral due to the self noise terms is shown to be inde- 
pendent of the mean flow. This contribution, as shown 
by Ribner [6], is primarily due to terms like ViVjV^v^. 
As explained by Mani [7], the distinction between shear 
noise and self noise becomes less important when the 
interpretation of the source terms is based upon Lil- 
ley’s [8] formulation rather than Lighthill’s equation. 
Assuming a normal joint probability for turbulent ve- 
locity components, as suggested by Ribner [6] and us- 
ing the isotropic turbulence model of Batchelor [9], the 
fourth-order correlation can be written as a linear com- 
bination of second-order correlations. 

ViVj v' k v[ - 

( 3 ) 

Further b y ass uming that the two-point velocity cor- 
relations, are separable in space/time factors, as 
postulated by Ribner, 

Viv'j = Rij(^)g(r) (4) 

the integration on separation vector of correlation can 
be carried out in closed form. The space factor, accord- 
ing to Batchelor, is written as: 

*.;(*) = T[(f + \u')f>a - \f'U,K\ ( 5 ) 


1. Source Spectrum Model 


The mean-square sound pressure in the far field, 
in absence of convection and refraction is written using 
LighthilFs acoustic analogy approach [4,5]. 


pW^) = 


Ri Rj Rk Rt 
16* 2 C^R 6 


U^ Ti “ 


The source strength is assumed to he dominated by 
unsteady momentum flux i.e. ~ pViVj. Vector £ 
in the above integration is the separation vector of the 
correlation between pV^Vj computed at y and 
at y + £ and r is the time delay of correlation. The 


where T = 3WW is the intensity of turbulence replac- 
ing the axial turbulence intensity in the Reichardt’s 
model and /' = df jd£. The three components of 
the separation vector of correlation are represented as 
( i for i = 1,2,3. Function /(£) can be expressed as 
/(£) = exp(-7r$ 2 /L 2 ). This expression makes /(£) de- 
crease to zero for large £ with sufficient rapidity to make 
Jo°° converge for m > 0. Length L is the 

longitudinal macroscale of turbulence, representing the 
linear extent of the region within which velocities are 
appreciably correlated. The time factor of correlation 
is expressed as g(r) = exp(—r/r c ) where r 0 denotes the 
characteristic time-delay in a moving reference frame. 
For axisymmetric jets, Davies et al. [10] have shown 
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that t c is proportional to the inverse of mean shear i.e. 
\Jt 0 — (< dU/dr ) where U denotes the mean velocity 
and r the radial distance. In addition, turbulence in- 
tensity can be expressed in terms of L and r 0 according 
to r 0 ~ L/y/T. Since the eddy length scale L t is related 
to kinetic energy of turbulence k = and its dis- 

sipation rate c as L t = l fc 3 ^ 2 /e, assuming that L e and L 
are proportional, it can be concluded that 1 jr 0 ~ c/k. 
The above conclusion can also be derived from the alge- 
braic turbulence model. Comparison of r 0 , as computed 
from k and e, with the corresponding results obtained 
from the velocity gradient and Reichardt’s aerodynamic 
solution are presented under the numerical discussions. 

Ribner [6] evaluates the contribution to self noise 
for various source strength components. The corre- 
sponding spectrum is evaluated by using a Fourier 
transform on the time delay of the correlation [2]. In 
the final analysis the source intensity spectrum can be 
expressed as 

J(fi)~ p 2 k^(Qr 0 ) i e'^? (6) 

To compute the noise field, this equation will later be 
used in conjunction with the refraction effect of the 
mean flow. The doppler effect relating the source fre- 
quency Q and the observed frequency / is expressed 
through a convection factor, C, as Q = 2 irfC. The 
singular behavior of the eddy convection factor C — 
(1 - M c cos6) at supersonic convection speeds (M c > 1) 
has been discussed by Ffowcs Williams [11]. The auto- 
correlation function in this region, which corresponds 
to emission of eddy Mach wave, becomes zero and thus 
the integrand remains finite. In our analysis a modified 
convection factor described by 

C = y/(l - M c cos0) 2 + (a c Jk VCoo) 2 (7) 

has been selected for numerical computation. The em- 
pirical convection constant a c was determined from 
comparison of prediction with data. A value of a c = .5 
is generally a good approximation for most observation 
angles. The convection Mach number was expressed as 
a function of the weighted average of the nozzle exit 
Mach number and local Mach M 

M c = .5 M + PcMj (8) 

Values in the range of .25 to .3 for the convection con- 
stant j3 c appear to yield best results. 

2, Sound/Flow Interaction 

Lighthill’s acoustic analogy approach, based on 
classical wave equation doesn’t incorporate the effect of 


surrounding mean flow on the sound radiation caused 
by convected multipole sources. It is clear that these 
pressure fluctuations should propagate through a region 
of nonuniform velocity and temperature before reaching 
the observer point. Thus the location of source within 
the jet plume has a strong effect on the amount of rar- 
diated sound. The mean flow results in not only the 
refraction of the radiated sound, but also an additional 
convection amplification due to fluid motion. A number 
of investigators have studied the radiation field of mul- 
tipole sources immersed in parallel sheared flows. Mani 
[7,12,13] studied the mean flow interaction of round jets 
for slug flow profiles assuming quadrupole sources con- 
vening along the centerline of the jet. Gliebe and Balsa 
[14] and Goldstein [15] extended the above analysis to 
arbitrary velocity profiles in high frequency and low 
frequency limits respectively. The generalization to ar- 
bitrarily located sources in continuously varying mono- 
tonic profiles were derived by Balsa [16] and Goldstein 
[17,18]. For supersonic jets, the high frequency solu- 
tion provides an adequate approximation as suggested 
by Tester and Morfey [19] and Pao [20]. 

Lilley’s equation is considered as the starting point 
by most recent investigators studying the sound/flow 
interactions. For a parallel flow, the Green’s function 
solution to a convected monopole of frequency Q is ob- 
tained by applying a sequence of Fourier transforms 
tions. The solution corresponding to multipole singu- 
larities are obtained by differentiating the monopole 
solution with respect to appropriate source coordi- 
nate; once in the case of a dipole source and twice 
for quadrupole sources. Balsa and Mani [2] provide 
a comprehensive analysis of the shielding effect of an 
axisymmetric mean flow, i.e. when the jet velocity and 
temperature profiles are arbitrary functions of the ra- 
dial variable. To study the effect of asymmetry in the 
mean flow, Goldstein [21] solves Lilley’s equation for 
high frequency multipole sources in a jet flow whose 
Mach number and temperature are functions of cross- 
flow coordinates. This leads to an ordinary differen- 
tial equation that traces the cross-stream projection of 
an acoustic ray starting from any given source loca- 
tion. The differential equation is solved repeatedly as 
one selects initial departure angle for rays leaving the 
source. Subsequent ray tracing through the shear layer 
and study of the spreading or focusing of the acous- 
tic emission as it approaches far field, results in a cir- 
cumferential directivity factor for each source volume 
element. 

For our axisymmetric geometry, Mani and Balsa’s 
formulation, as used in the MGB code, has been applied 
whereby the turbulent properties of the jet has been 
coupled with its acoustic radiation. Thus the mean 
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square pressure in the far field is expressed as 

p 2 (R } 0 } Q) = I A (a xx + 4a xy -j- 2 a yy + 2 a yz )dy (9) 
Jg 

The contribution to noise field due to each of the 
quadruples contained within a turbulent eddy vol- 
ume element is described by the directivity factors 
a xr ,...,a yz . The corresponding weighting factors were 
derived by Ribner [6]. A is a factor related to the source 
intensity and frequency 

j ^ ^ \ p / \ t /fo) 

~ {AifRC^Cfil- Mcosef{\- M c cos6f K ’ 


The directivity factors are functions of the shielding 
function g 2 which is expressed in terms of the polar 
observation angle 0, local sound velocity C and Mach 
numbers M(r) = U{r)/Coo and M c = U c /Coo ■ 


2 _ ( 1 - AQ 2 (%0 -cos?0 

M"; — 2 


(1 — M c cos0 ) 


(ii) 


The mean flow variables, as obtained from CFD, are 
used in computation of g 2 (r) and hence the directivity 
factors a ir ,.,.,a yi . Location r where g 2 (r) changes in 
sign is called transition point or turning point. When 
g 2 (r) is negative between the source and observer, the 
possibility of fluid shielding exists. The position of the 
source with respect to the turning points of g 2 {r) con- 
tributes to the amount o f fluid shielding through factor 
exp(— 2K \J\ g 2 (r) \ dr) . The limits of the integra- 
tion are determined by the location of the source rela- 
tive to the turning points. 

Wave number K = u>/C in the above exponen- 
tial factor tends to increase the attenuation of the high 
frequencies. An extensive list of directivity factors cor- 
responding to a case with two turning points have been 
listed in [2]. A maximum of four turning points have 
been considered in the present computations. A correc- 
tion of one Doppler factor has been utilized for source 
volume effects according to Ffowcs Williams [11]. The 
correction for flight speed has been included using flight 
dynamic factor (1 4* M^cosO)" 1 where M©© is the flight 
Mach number. 


Numerical Results and Discussion 

To obtain the flow field solution a 141 x 60 grid 
(Figure 1), highly clustered around the jet lip-line 
was used. A run ti me of over four hours on CIRAY- 
YMP was required to achieve the necessary conver- 
gence within PARC-ke code. The aerodynamic results 
obtained from PARC code were input to a noise pre- 
diction code. Aerodynamic and acoustic calculations 


were compared with available data. Also turbulent 
intensity and characteristic time delay of correlation, 
as predicted by PARC code, were compared with Re- 
ichardt’s solution. All laser velocimeter, shadowgraph 
and acoustic tests for this model were conducted in the 
General Electric anechoic free-jet facility [1]. The test 
point (corresponding to the design condition) is listed j 

in table 1. The Mach contour plot (Figure 2) indicates ! 

the for mati on of a Mach disc inside the nozzle as well 
as a weak shock outside the exit plane, also seen in the J 

sh adowgraph photo [2] . 

Figures 3 and 4 demonstrate the close agreement \ 

with data in the velocity profiles on the lip-line, cen- * 

terline, and also at four axial positions along the jet. ! 

The small oscillations in velocity near the centerline 
can better be predicted by increasing the grid resolu- 
tion from exit to 10 diameters downstream. Figure 5 I 

exhibits the contour plot for the intensity of turbulence, : 

indicating the extension of noise sources to nearly 18 di- | 

ameters downs tre am. The contour levels are nondimen- 
sionalized with respect to reference sound speed, a 2 r , 1 

which is computed at the stagnation temperature. Fur- § 

ther details of the turbulence intensity distribution are f 

demonstrated in Figures 6 and 7. The axial distribu- | 

tion of the turbulence intensity near the lip-line (Figure 
6a) shows a level of 13 percent near the exit plane with 
a growing decay rate after 12 diameters, generally in 
good agreement with data. Reichardt’s model predicts 
a much faster decay after 5 diameters. The radial dis- 
tribution (Figure 6b) was plotted only at X/D = 8.21 
where' test data was available for comparison. Obvi- 
ously PARC code does a better job in predicting the 
maximum turbulence level. Near the centerline, which 
is still within the potential core, the CFD computations 
clearly under-predict the data by 3 to 4 percent. This 
can be attributed to the boundary conditions used in 
a CFD solution and should not introduce any signif- 
icant error in noise computation because of low level 
of source intensity and high shielding effect. Further 
comparison of the source term using the two aerody- 
namic prediction methods is illustrated in Figure 7. 
ReichardVs model predicts a faster spreading rate and 
also fails to properly predict the ending of the potential 
core where increased turbulence intensity near the cen- 
terline is expected. In addition, Reichardt’s model does 
not accurately predict the axial decay of the turbulence 
intensity (see Figure 6a). 

The characteristic time delay r 0 was obtained us- 
ing two ^different definitions alluded to earlier. The cor- 
responding Strouhal numbers defined as ( t/k)(D/Uj ) 
an(l (dU/dr)(D/Uj) are plotted in Figures 8 and 9. 

It is clear that outside the proportionality factor ay 
needed for the computation of r C) these curves demon- 
strate remarkable resemblance, showing a peak value 
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moving in toward the centerline with increasing X. A 
third curve was also included based on Reichardt’s so- 
lution, showing the same general features but failing to 
properly predict the radial location of the maximum. 

In our acoustic computations r Q was calculated 
from 1 /t 0 = aje/k. Factor aj can shape the sound 
spectrum in terms of the location of the maximum 
and should be determined empirically. A value of 
aj = 2 was selected in these calculations. For non- 
axisymmetric geometries, additional modification in 
the definition of t 0 is necessary as the circumferential 
gradient of the mean velocity also contributes to this 
factor. The dissipation rate of turbulent kinetic energy 
is plotted in Figure 10. The contour levels are nondi- 
mensionalized with respect to a 3 r /L r where L r is the 
reference length used in the PARC code. The acoustic 
computations for a 40 foot radius are presented in Fig- 
ures 11 and 12. The overall sound pressure (OASPL) 
directivity as obtained from PARC shows reasonable 
agreement with data at most observation angles. The 
empirical factors /3 C and aj were selected to obtain the 
best fit at angles where the jet mixing noise is max- 
imum. Obviously for off-design conditions, the shock 
noise dominates at angles close to upstream axis and 
several dB’s difference in our predictions would not sig- 
nificantly alter the directivity of the total noise. The 
over-prediction of up 5 dB’s at 6 = 80 degrees can 
be improved by selecting an angle dependent convec- 
tion constant a c . Reichardt’s solution also generates 
results in close agreement with data as the empirical 
constants used in this aerodynamic prediction model 
are calibrated for a conical flow field. The spectral 
components of noise, based on a one-third octave band 
are shown in figure 12. The overall agreement seems 
acceptable. 


Concluding Remarks 

The source strength has been successfully pre- 
dicted using PARC code with a k-e turbulence model. 
The resulting noise predictions compare reasonably 
with experimental data. This application has removed 
four empirical constants used in the Reichardt’s model, 
which are calibrated for a conical flow field. The limita- 
tions on aerodynamic grid selection has been removed 
by adopting a two-stage aerodynamic and acoustic al- 
gorithm. Further, the time-delay of correlation has 
been computed using two different models. The acous- 
tic computations, according to the present model, is 
based on two fundamental assumptions as suggested 
by Ribner and Batchelor. A proper prediction of the 
noise spectrum depends on an accurate assessment of 
the time-delay of correlation. 

For non-axisymmetric geometries, a proper cir- 
cumferential averaging of the flow variables can be used 


as a first step in the computation of the polar directiv- 
ity of jet noise. Also contributions from other compo- 
nents of the mean velocity gradient should be included 
in the computation of t 0 . More elaborate analysis such 
as geometric ray theory is needed to predict the az- 
imuthal spreading of jet noise. In addition, the empiri- 
cal constants introduced in the determination of super- 
sonic convection factor should further be investigated 
for other geometries. 
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Nomenclature 

local and ambient speeds of sound 
convection factor 
exit and throat diameters 
observed frequency 
source intensity spectrum 
turbulent kinetic energy = vJvT/2 
wave number 

longitudinal macroscale of turbulence 

eddy length scale 

local and exit Mach numbers 

convection and flight Mach numbers 

radial distance 

source to observer distance 

observer coordinate 

intensity of turbulence = t£vJ/3 

unsteady momentum flux 

mean velocity in direction of flow 

exit jet velocity 

fluctuating velocity component in X* 
direction 

local velocity component in X, 
direction 

two-points time-delayed second-order 

correlation = u,(y, t)vj (y + £ t + r)dt 

source coordinate 

convection constants 

proportionality factor (1 /r 0 = a/e/fc) 

Kronecker’s delta 

dis sipation rate of kineti c energy 

= i '(dvifdXjXdvi/dXj) 

polar angle with respect to 

downstream axis 

source factor 

vector separation of correlation 
local and ambient densities 
molecular viscosity 
time delay of correlation 
characteristic time delay of correlation 
azimuthal angle 

observer and source radian frequencies 
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Table 1. Design Parmaters for C-D Nozzle 

Throat Diameter 
Exit Diameter 

Distance from Throat to Exit 
Exit Velocity 
Ambient Velocity 
Pressure Ratio 
Stagnation Temp. 


5.1 in. 
5.395 in. 
5.525 in. 
2409 fps 
400 fps 
3.121 
1716*R 
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Fig. 2: Mach number contour plot ( Mj = 1.4). 
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Fig. 7: Comparison of PARC-ke turbulent intensity profiles with Reichardt’s 
solution at different axial locations, (a) PARC-ke (b) Reichardt’s model. 
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Fig. 3: Strouhal No., (a) St. defined as (e/k)(D /Uj) 
obtained from PARC-ke. (b) St. defined as 
( dU/dr)(D/Uj ) obtained from PARC-ke , (c) 
same as (b) but obtained from Reichardt’s 


Fig. 9: Strouhal No., (a) and (b) obtained from 
PARC-ke, (c) Reichardt. 










Fig. 10: Contour plot for dissipation rate of turbulent kinetic energy (PARC-ke). 
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Fig. 11: Comparison of the overall sound pressure level directivity with data [1] 
on a 40 foot radius. 
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Fig. 12: Comparison of the spectral components of noise with data [1] at different 
polar angles based on a 1/3 octave center frequency (band number 11 
corresponds to 50 Hz.). n 
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Fig. 12: Comparison of the spectral components of noise with data [1]. 
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